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ABSTRACT:  Numerical  simulators  of  ground- water  flow  and  transport  are  frequently  used 
to  determine  streamlines  and  to  estimate  travel  times  and  pathways  of  contaminant  move¬ 
ment.  These  items  are  often  obtained  by  tracking  conceptual  water  particles  through  the 
computational  grid  using  model-calculated  flow  velocities.  Such  tracking  is  also  an  important 
component  of  the  Lagrangian  part  of  many  methods  for  transport  modeling,  in  particular 
the  Eulerian-Lagrangian  localized  adjoint  method  (ELLAM).  A  procedure  for  exact  analyt¬ 
ical  particle  tracking  is  presented,  given  a  lowest-order  Raviart-Thomas  velocity  field  v  on 
a  rectangular  spatial  grid,  with  linear  temporal  interpolation  of  v  from  the  beginning  to 
the  end  of  a  time  step.  This  includes  xt-  and  yt-bilinearity  in  the  x-  and  y-components, 
respectively,  of  v.  Previous  authors  assumed  that  v  was  steady,  or  that  its  time  derivative 
was  constant  in  space.  Transience  in  v  allows  a  particle  to  reverse  its  direction  during  a 
time  step.  The  added  effect  of  bilinearity  can  be  significant,  especially  when  v  varies  in 
time  due  to  changes  in  well  pumping  rates  or  variable  recharge.  These  effects  are  discussed 
qualitatively  and  illustrated  with  test  problems  that  compare  the  accuracy  of  the  tracking 
methods. 


1  INTRODUCTION 

Problems  of  ground-water  flow  and  solute  transport  have  been  intensely  studied  in  recent 
years  due  to  increased  awareness  of  the  susceptibility  of  ground  water  to  contamination.  Nu¬ 
merical  simulators  of  ground-water  flow  and  solute  transport  are  common  tools  employed  in 
these  analyses.  Numerical  models  are  frequently  used  to  determine  streamlines  and  estimate 
travel  times  and  pathways  of  contaminant  movement.  These  items  are  often  obtained  by 
tracking  imaginary  particles  of  water  through  the  model’s  computational  grid  using  model- 
calculated  flow  velocities.  Pollock  (1988,  1989)  used  particle  tracking  to  determine  the  loca¬ 
tion  of  streamlines  and  to  estimate  the  time  required  for  water  to  traverse  different  segments 
of  a  model  domain.  Tracking  was  also  applied  to  estimate  the  extent  of  advective  transport  of 
solutes  (Garabedian  &  Konikow  1983).  Other  uses  include  its  incorporation  into  numerical 
models  of  solute  transport  that  are  based  on  the  advection-dispersion  equation  (Konikow 
&  Bredehoeft  1988;  Prickett  et  al.  1981;  Zheng  1989;  Healy  &  Russell  1993).  This  paper 
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presents  a  new  method  for  particle  tracking. 

The  first  step  in  a  general  particle  tracking  procedure  is  to  solve  the  ground-water  flow 
equation  using  standard  finite-difference  or  finite-element  space  and  time  grids  to  obtain 
estimates  of  hydraulic  head  or  pressure  at  fixed  node  locations  and  times.  Darcy’s  equa¬ 
tion  is  then  employed  to  calculate  velocities  from  the  head  or  pressure  field.  Alternatively, 
heads  and  velocities  can  be  solved  simultaneously  as  in  the  method  of  mixed  finite  elements 
(Raviart  &;  Thomas  1977;  Russell  &  Wheeler  1983).  Either  approach  results  in  velocities 
that  are  calculated  at  the  interfaces  between  adjacent  nodes.  These  velocities  are  then  used 
to  calculate  the  pathlines  and  travel  times  of  particles. 

We  assume  a  Raviart-Thomas  (1977)  velocity  field:  within  each  spatial  cell,  the  x- 
component  is  continuous  piecewise-linear  in  the  ^-direction  and  discontinuous  piecewise- 
constant  in  the  ^-direction,  and  the  y-component  is  the  reverse.  This  was  also  assumed  by 
Pollock  (1988,  1989)  and  Schafer- Perini  &  Wilson  (1991).  This  vector  field  lends  itself  nat¬ 
urally  to  linear  interpolation  for  determining  velocities  in  the  interior  of  each  cell.  Linear 
interpolation,  however,  does  produce  discontinuities  in  the  velocity  field  at  cell  interfaces 
(Goode  1990). 

Calculation  of  travel  path  lines  for  any  particle  requires  integration  of  the  particle  velocity 
over  time.  This  integration  can  be  accomplished  by  one  of  three  approaches:  analytical, 
numerical,  or  semianalytical.  The  analytical  approach  (Javandel  et  al.  1984)  produces  exact 
solutions,  but  only  for  a  limited  number  of  ideal  cases  of  steady  flow,  homogeneous  media, 
and  simple  geometry.  Such  cases  seldom  arise  in  practice. 

Numerical  techniques  for  integration  include  the  explicit  single-step  method  (Goode 
1990),  the  first-order  Euler  (Lu  1994),  and  the  fourth-order  Runge-Kutta  (Nelson  1978; 
Shafer  1987;  Zheng  1989).  These  schemes  are  not  limited  by  transient  velocities  or  complex¬ 
ity  in  the  velocity  fields.  The  single-step  explicit  method  is  computationally  simple,  but  of 
limited  accuracy.  The  other  two  methods  can  attain  a  high  degree  of  accuracy,  but  may 
require  a  large  number  of  steps  (and  therefore  computation  time)  to  do  so. 

The  semianalytical  technique  combines  aspects  of  analytical  and  numerical  methods.  It 
makes  use  of  an  analytical  solution  to  the  integral  within  an  individual  space  and  time 
cell  under  the  assumption  of  the  Raviart-Thomas  velocity  field.  Tracking  is  then  conducted 
through  one  cell  at  a  time.  This  idea  was  first  presented  by  Pollock  (1988,  1989)  for  steady- 
state  flow  conditions  and  has  since  been  used  by  Schafer-Perini  &  Wilson  (1991)  and  Healy 
&  Russell  (1993).  It  was  extended  to  include  nonsteady-state  conditions  under  the  added 
assumption  that  velocities  varied  linearly  with  time  within  each  time  step  and  that  the 
spatial  derivative  of  velocity  was  constant  in  time  (i.e.,  d2v/dxdt  =  0)  (Lu  1994). 

In  this  paper  we  extend  the  work  of  Lu  (1994)  and  obtain  the  tracking  equations  for 
nonzero  d2v/dxdt.  Test  problems  are  presented  to  demonstrate  the  sensitivity  of  results  to 
inclusion  of  this  term.  Because  of  the  complex  nature  of  the  tracking  equations,  a  detailed 
algorithm  is  provided.  The  algorithm  permits  straightforward  treatment  of  single  or  double 
reversals  in  velocity  direction  that  may  occur  within  a  single  time  step. 

2  TRACKING  EQUATIONS 

2.1  Decoupled  ordinary  differential  equations 

Let  C  —  [xi- 1,  Xi]  x  [y-j-i- yj]  x  [tn.  tn+1]  be  a  space-time  cell  in  a  two-dimensional  rectangular 
grid  over  a  time  step.  The  developments  described  below  will  generalize  in  an  obvious  way 
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to  three  space  dimensions.  In  C,  suppose  that  a  particle  is  located  at  (xo7  ijq)  at  time  to.  Let 

v0  =  «,  vyo )  =  (■ vx(x0 ,  y0,  to),  ^(xo,  y0,  *o))  (1) 

be  the  velocity  vector  at  this  location  and  time.  Spatially,  the  velocity  belongs  to  the  lowest- 
order  Raviart-Thomas  space  RT(h  and  the  RT0  coefficients  are  linear  functions  of  time,  so 
that  vx  is  bilinear  in  x  and  t  and  constant  in  y.  and  vy  is  bilinear  in  y  and  t  and  constant  in 
x.  Accordingly,  if  we  let 


Vx  = 

dvX  (t  l 
dx  {to) 

(independent  of  x  and  y), 

(2) 

Vy  - 

dvv  (t  1 
dv  (k) 

(independent  of  x  and  y), 

(3) 

^  = 

dvx 

dt{Xo) 

(independent  of  y  and  t). 

(4) 

vyt  - 

dvy  .  . 

at  w 

(independent  of  x  and  t), 

(5) 

Vxt  — 

d2vx 

dxdt 

(constant), 

(6) 

vyt  ~ 

d2vy 

dydt 

(constant), 

(7) 

then  we  can  write  the  velocity  components  on  C  as 

Vx(x ,  t )  =  v%  +  vxAt  +  ( vx  +  vxtAt)Ax7  (8) 

vv(y,t)  =  v%  +  v%At+(vy  +  vytAt)Ay,  (9) 

where  Ax  —  x  —  x0  and  Ay  —  y  —  y0. 

The  tracking  of  the  particle  across  C ,  which  amounts  to  determining  the  trajectory 
(x{t).  y(t)),  is  governed  by  the  ordinary  differential  equations 

x'{t)  -  vx(x(t),t)7  x(t0)  -  x o,  (10) 

y'(t)  -  vv(y(t),t),  y(to)=y0.  (11) 

Because  vx  does  not  depend  on  y  and  vy  does  not  depend  on  x,  these  ODE’s  are  decoupled, 
and  each  can  be  solved  on  C  without  reference  to  the  other.  The  analogous  observation  holds 
in  three  dimensions. 

2.2  Global  tracking  algorithm 

The  global  tracking  algorithm  consists  of  a  sequence  of  local  tracking  steps,  each  of  which 

is  confined  to  a  single  space-time  cell  C.  A  step  begins  at  an  initial  time  t0  >  tn,  which  is 
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either  the  starting  time  for  the  trajectory  or  the  time  at  which  the  trajectory  enters  C.  The 
step  lasts  until  a  final  time  tf  <  tn+1,  which  is  either  the  ending  time  for  the  trajectory  or 
the  time  at  which  it  leaves  C  (note  that  it  leaves  at  time  tn+1  if  it  does  not  cross  the  spatial 
boundary  of  C ). 

Because  the  ODE’s  (10)  and  (11)  are  decoupled,  one  can  separately  determine  times 
tx  >  t0  and  ty  >  t0  at  which  x(t)  and  y(t).  respectively,  leave  the  ranges  corresponding  to  C. 
Then  tf  is  set  equal  to  the  minimum,  tf  —  min {4,  ty}.  Similarly,  in  three  dimensions,  tf  is 
the  minimum  of  three  times  determined  by  independent  one-dimensional  calculations.  Thus, 
the  objective  of  the  detailed  algorithm  below  is  to  carry  out  one  local  step  in  one  dimension. 
These  local  one-dimensional  steps  can  be  compared  to  find  f/,  then  concatenated  to  follow 
the  global  trajectory. 


2.3  Global  tracking  equations 

Accordingly,  we  confine  ourselves  henceforth  to  a  single  dimension  x.  The  ODE  (10),  with 
terms  defined  in  (2),  (4),  (6),  and  (8),  can  be  solved  analytically.  We  drop  the  superscript  x 
from  Vq  and  vf  in  (1),  (4),  and  (8).  Suppose  for  now  that  these  definitions  hold  globally,  i.e., 
ignore  the  fact  that  they  are  different  on  cells  other  than  C.  The  analytical  solution  is 


x(t)  =  xq  -\ — -  (e2 

vxt  v 


Vt  ri'.,  At 


-i) 


+  U’o  - 


Vxt  J 


)e^{M+z)2  rz 


2  V. 


X 


xt 


X 


erfc  ( \l -Ai_  )  _  erfc  ^  A t  +  ^ 
2  vxt 


2  V  vxf 

Vxt  >  0, 


X 


(t)  =  Xo  +  —  UVxiAt2+VxAt  -  l) 
Vxt 


(12) 


+  A)  - 


vxvt 


/  V  -Vxt 


X 


X 


/ 


Vxt 


vxt  Vx 
2  vxt 


es  ds ,  vxt  <  0, 


,  ,  eVxAt  -  1  eVxAt  -  (1  +  vxA t) 

x(t)  =  X0  +  Vo - h  Vt - 2 - : 


(13) 


Vxt  =  0,  vx  ±  0, 

x{t)  =  x0  +  v0A  t  +  ^vtA  t2,  vxt  =  vx  =  0. 


(14) 

(15) 
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Equations  (14)  and  (15)  are  those  of  Lu  (1994).  In  our  more  general  context,  (12)  through 
(15)  are  not  satisfactory  for  a  practical  implementation,  due  to  difficulties  with  roundoff 
error.  Small  values  of  vxt  require  an  alternative  equation  that  eases  the  transition  from  (12) 
or  (13)  to  (14).  Let  e  denote  the  machine  precision,  i.e.,  the  smallest  positive  number  such 
that  1  +  e  >  1.  Analysis  of  (12)  and  (13)  finds  roundoff  error  0(v~le).  Then  asymptotic 
techniques  to  determine  the  behavior  of  x(t)  for  small  vxt  lead  to  (14),  plus  a  first-order  term 
in  vxt ,  with  truncation  error  0(u2t).  As  vxt  tends  to  0,  the  roundoff  error  increases  and  the 
truncation  error  decreases,  yielding  a  threshold  \vxt\  =  6  above  which  (12)  or  (13)  is  used, 
and  below  which  (14)  with  the  additional  first-order  term  is  used.  The  threshold  is  where 
the  two  errors  are  equal,  0{v~fc)  —  0(vxt),  so  that  6  —  0(e1//3)  and  the  error  is  0(e2//3).  The 
actual  implemented  equations  are 


x{t )  =  right-hand  side  of  (12),  vxt  >  <5, 

x{t )  =  right-hand  side  of  (13),  vxt  <  — 5 , 

,  ,  eVxAt  -  1  eVxAt  -  (1  +  vxA t) 

x(t)  =  Xq  +  Vo - 1-  vt - 5 - 


(16) 

(17) 


+  Vxt 


V°  t1  At2eVxM  -  eVxM  ~  0- +  v*AtY 


vTj\2 


where 


-3 


ev*At  -  (1  +  vxAt  +  \{vxA tf  +  At)3)' 


\vxt\  <  5,  vx  not  small, 


5  = 


1/3 


X 


2evt  max{l,  e'JxAt} 

vtVx At4  e^At  -  (1  +  v^A t  +  |K At)2) 


■( 


+ 


u0At4  eVxAt  -  1 

_L 

VfA  t6 

8  vx 

48 

-1/3 


(18) 


(19) 
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In  (18)  and  (19),  if  vx  is  small,  make  the  following  replacements,  letting  s  —  vxA t: 

Replace  with  if  |s|  < 

Ai(l  +  f )  (he)1/* 

At2(±  +  §)  (24c)1/4 

At3(i  +  |)  (f  J)1'5 

A*4(A  +  5o)  <720e)1/6 
At3(l  +  £)  (120c)1/5 

The  next  section  discusses  how  these  equations  are  used  within  the  confines  of  a  space-time 

cell  C. 

3  LOCAL  TRACKING  ALGORITHM 

A  complication  in  determining  the  time  tf  when  a  trajectory  leaves  a  cell  is  the  possibility 
of  reversals  in  the  direction  of  the  trajectory.  Thus,  one  cannot  simply  assume  that  if  the 
trajectory  is  within  the  range  of  C  at  the  initial  and  final  times,  then  it  must  be  within  the 
range  at  all  intermediate  times.  This  section  analyzes  the  possibilities  for  reversals  and  then 
bases  a  detailed  algorithm  on  the  results. 

3.1  Reversals  of  direction 


Equation  (8)  for  the  velocity  v(x,  t)  =  vx(x,  t )  can  be  rewritten  in  the  form 

v(x, t )  =  vxt  (^Ax  +  (a t  + 

Vxt' 

.  (  vxvt\ 

+  V(l 

Vxt  7^  0, 

(20) 

\  Vxt  ) 

v(x,  t )  =  vxAx  +  vtA  t  +  u0, 

c2 

a 

C-+. 

II 

CD 

(21) 

where  x  =  xo  +  Ax,  t  =  to  +  At,  and  we  think  of  Ax  and  At  as  the  variables.  Reversals 
can  occur  if  v  =  0  somewhere  in  C,  so  we  concern  ourselves  with  the  locus  of  points  in  the 
(Ax,  At)-plane  where  v  vanishes. 

If  vxt  —  0,  then  by  (21)  this  locus  is  a  straight  line  in  the  plane.  If  this  line  meets  C 
and  has  negative  slope  (vx  and  vt  have  the  same  sign),  the  two  essentially  distinct  reversal 
situations  are  depicted  in  Fig.  1.  Each  box  represents  the  space-time  cell  C  in  the  xt-plane  (x 
on  the  horizontal  axis),  the  solid  line  is  the  locus  of  v  —  0,  and  the  dashed  line  is  a  possible 
trajectory  with  slope  dt/dx  —  1/v  at  each  point.  In  the  left  figure,  vx  <  0  and  vt  <  0,  so 
that  v  >  0  in  the  lower  left  and  v  <  0  in  the  upper  right.  The  reverse  is  true  in  the  right 
figure.  Note  that  the  trajectory  can  only  cross  the  solid  line  vertically  (v  =  0  implies  that 
the  slope  1/v  is  infinite).  If  the  locus  has  positive  slope,  mirror  images  of  Fig.  1,  with  the 
signs  of  v  reversed,  cover  the  possibilities.  In  all  cases,  there  is  at  most  one  reversal,  and 
its  presence  can  be  detected  by  the  initial  and  final  velocities  v0  and  v f  having  opposite 
signs.  These  possibilities  arise  in  Lu  (1994).  As  a  special  case  of  (21),  if  the  velocity  is  steady 
(vt  =  0),  then  the  solid  line  in  Fig.  1  would  be  vertical,  and  the  trajectory  could  not  cross  it; 


ea— 1 

V,X  \ 
S~(l+s) 


i(iAtV-5tg±2>) 


~(l  +  g+|g2+^g0 
V% 

es-(l +s+^s2) 
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Figure  2:  Reversal  possibilities  with  bilinear  velocity  function,  v  —  0  on  intersecting  lines. 

in  other  words,  reversals  cannot  occur  in  a  steady  velocity  field,  as  in  Pollock  (1988,  1989) 
and  Schafer-Perini  &  Wilson  (1991). 

Next,  suppose  that  vxt  ^  0,  so  that  (20)  applies.  The  locus  of  v  =  0  is  a  hyperbola 
with  vertical  asymptote  Ax  =  —vt/vxt  and  horizontal  asymptote  At  =  —vx/vxt.  unless 
Vo vxt  —  vxvt.  in  which  case  the  locus  consists  of  the  two  asymptotes.  The  latter  case  is 
depicted  in  Fig.  2,  with  vxt  <  0  in  the  left  figure  and  vxt  >  0  in  the  right.  In  either  situation, 
at  most  one  reversal  is  possible,  and  the  signs  of  the  initial  and  final  velocities  v0  and  v /  will 
flag  it. 

Now  suppose  that  v0vxt  <  vxvt.  so  that  the  locus  is  a  hyperbola  whose  branches  lie  in 
the  first  and  third  quadrants  relative  to  the  asymptotes  and  have  negative  slope.  The  locus 
divides  the  plane  into  three  regions.  The  possibilities  can  be  classified  according  to  the  sign 
of  vxt,  which  determines  the  sign  of  v  in  each  region,  and  according  to  which  regions  the 
lower-left  and  upper-right  corners  of  C  lie  in.  The  six  essentially  distinct  cases  are  shown  in 
Fig.  3,  with  vxt  <  0  in  the  three  left  figures  and  vxt  >  0  on  the  right.  In  all  but  one  of  the 
cases,  at  most  one  reversal  is  possible,  and  it  can  be  detected  as  above.  The  possibility  of 
two  reversals  arises  in  the  middle-left  picture  in  Fig.  3;  however,  it  can  be  flagged  by  the  fact 
that  the  initial  and  final  velocities  are  negative  (vq  <  0,  v /  <  0),  while  the  net  movement 
from  the  initial  to  the  final  position  is  in  the  positive  (opposite)  direction  (xf  >  xq).  The 
cases  for  vovxt  >  vxvt  are  covered  by  mirror  images  of  Fig.  3  with  the  signs  of  v  reversed. 

The  reversal  cases  can  be  summarized  as  follows.  At  most  two  reversals  are  possible.  If 
there  are  two,  then  Xf  —  x0  has  sign  opposite  to  that  of  v0  and  Vf.  and  conversely.  If  there 
is  one,  then  v0  and  Vf  have  opposite  signs,  and  conversely. 
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Figure  3:  Reversal  possibilities  with  bilinear  velocity  function,  v  =  0  on  hyperbola. 

3.2  Detailed  local  tracking  algorithm 

The  algorithm  described  here  is  based  on  the  analytical  equations  and  the  reversal  pos¬ 
sibilities  derived  above.  The  objective  is  to  detect  the  earliest  time  at  which  a  trajectory 
leaves  the  space-time  cell  C.  If  there  is  no  reversal  during  the  time  step,  and  the  trajec¬ 
tory  crosses  the  cell  boundary,  this  leaving  time  is  found  by  solving  the  nonlinear  equation 
x(t)  =  boundary  coordinate,  where  x(t)  is  given  by  (16),  (17),  or  (18).  The  solution  is  ob¬ 
tained  by  Newton’s  method,  using  (8)  for  the  derivative  x'(t )  =  vx(x(t).  t),  with  a  bisection 
method  as  a  backup  if  a  Newton  step  fails  to  make  x(t)  closer  to  the  desired  value.  If  there 
is  a  reversal,  then  its  time  t*  is  found  by  solving  v(x(t).  t.)  —  0,  with  v  given  by  (8).  Again 
Newton/bisection  is  used,  with  derivative 

v'(t)  =vt  +  vxtAx  +  (vx  +  vxtA  t)v(t).  (22) 

Once  t*  is  found,  it  is  known  that  there  are  no  reversals  between  t0  and  t*.  so  that  the  time 
interval  up  to  t*  can  be  treated  as  a  step  without  reversals. 

In  the  following  step-by-step  description,  t0  is  the  current  time,  at  which  the  local  tracking 
step  begins;  x0  =  x(t0)  is  the  current  position,  from  (16)  through  (18);  Vq  is  the  current 
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velocity,  from  (1);  vx,  vt,  and  vxt-  from  (2),  (4),  and  (6),  are  the  current  parameters  in  the 
velocity  equation  (8),  and  may  depend  on  the  velocity  direction  if  xo  is  at  a  cell  boundary 
(choose  those  corresponding  to  the  cell  being  entered  by  the  trajectory);  tf  is  the  end  of  the 
current  time  step,  the  goal  of  the  local  tracking  step;  and  xL  and  xR  are  the  left  and  right 
endpoints  of  the  cell  containing  Xq  (the  cell  being  entered  by  the  trajectory  if  x0  is  at  a  cell 
boundary).  Note  that  the  “cell  being  entered”  is  unambiguous,  because  v  is  continuous  at 
cell  boundaries. 

At  the  start  of  each  local  tracking  step,  to,  vq,  vx,  vt,  vxt,  tf,  xl,  and  xR  are  set.  To 
complete  the  step  and  be  ready  for  the  next  step,  the  following  algorithm  determines  new 
values  for  these  parameters. 

1.  Compute  the  final  position  Xf  —  x(tf)  from  (16)  through  (18)  and  the  final  velocity 
Vf  =  v(xf,tf )  from  (8),  as  if  the  current  velocity  parameters  were  global. 

2.  If  vo  and  v /  have  opposite  signs,  go  to  step  6.  (If  vq  —  0,  assign  the  sign  of  vt\  if  v /  =  0, 

dv 

assign  the  sign  of  ~—(xf, tf)  =  -(■ vt  +  vxt(xf  -  x0 )).) 

3.  [uo  and  v j  have  the  same  sign]  If  vq  and  Vf  are  negative,  go  to  step  5. 

4.  [uo  and  Vf  are  positive]  If  xq  <  Xf  <  xR,  go  to  step  15.  If  xR  <  Xf,  set  t\  —  to,  £2  =  tf 

and  go  to  step  16.  If  Xf  =  xR,  go  to  step  18.  Otherwise  (xf  <  x0),  go  to  step  9. 

5.  [u0  and  Vf  are  negative]  If  xr  <  Xf  <  x0,  go  to  step  15.  If  Xf  <  xl,  set  t\  —  t0,  t2  —  tf 

and  go  to  step  17.  If  Xf  =  xl,  go  to  step  20.  Otherwise  (xq  <  Xf),  go  to  step  12. 

6.  [there  is  one  reversal]  Set  t\  —  to,  h  —  tf,  and  use  Newton/bisection  with  (22)  to  find 

the  unique  t*,  t\  <  t*  <  t2,  such  that  v(x(t*),t*)  =  0.  Evaluate  x*  =  x(t*). 

7.  If  x*  <  xl,  set  t\  —  to,  t2  —  t*,  and  go  to  step  17.  If  xr  <  x*,  set  ti  =  to,  t2  —  t*,  and 
go  to  step  16. 

8.  [reversal  occurs  within  current  cell]  If  Xf  <  xr,  set  t\  —  t*,  t2  —  tf,  and  go  to  step  17. 
If  xr  <  Xf,  set  t\  —  t*,  t2  —  tf,  and  go  to  step  16.  If  xl  <  Xf  <  xr,  go  to  step  15.  If 
Xf  =  xl,  go  to  step  20.  Otherwise  (xf  =  xr),  go  to  step  18. 

9.  [two  reversals,  Vo  and  v /  positive,  net  movement  to  left]  Set  t.\  —  to,  t2  —  tf,  and  use 
Newton/bisection  with  (8)  to  find  the  unique  t,t\  <t  <  t2,  such  that  x(t)  =  (x0+Xf)/ 2. 
Then  set  t\  =  t0,  t2  =  t,  and  use  Newton/bisection  with  (22)  to  find  the  unique  t*, 
t\  <  t*  <  t2,  with  v(x(t*),  t*)  —  0;  this  is  the  time  of  the  first  reversal.  Set  x*  —  x(t*). 

10.  If  xr  <  x*,  set  t\  =  to,  t2  =  t* ,  and  go  to  step  16. 

11.  [first  reversal  occurs  within  current  cell]  If  x[t)  <xl,  set  t\  —  t*,t2  —  t,  and  go  to  step 
17.  Otherwise  (x(t)  >  xL),  go  to  step  22. 

12.  [two  reversals,  v0  and  v /  negative,  net  movement  to  right]  Same  as  step  9. 

13.  If  x*  <  xl,  set  t\  —  to,  t2  —  t*,  and  go  to  step  17. 
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14.  [first  reversal  occurs  within  current  cell]  If  xr  <  x(t),  set  t\  —  t*,  t2  —  t,  and  go  to  step 
16.  Otherwise  (x(t)  <  xr),  go  to  step  22. 

15.  [tracking  for  full  time  step  stays  within  current  cell]  Set  to  —  tf,  tf  to  the  end  of 
the  next  time  step,  x0  =  x(t0),  v0  =  v(x0,t 0),  vx  =  —  (x0,t0),  vt  =  —  (x0,t£),  and 

cPv 

vxt  —  „  -  fa,  t|t),  where  tt  refers  to  the  next  time  step  that  begins  at  to;  xl  and  xr 
OxOt 

are  unchanged.  Go  to  step  23. 

16.  [tracking  crosses  right  endpoint  before  end  of  time  step]  Use  Newton/bisection  with 
(8)  to  find  the  unique  t' ,  t\  <  t'  <  t2,  such  that  x(t ')  =  xr.  Set  t0  =  t' .  x0  =  x(t')  =  xr , 

Ov  Ov  O^v 

v0  =  v(x0,t0),  Vx  =  —(x$,to),  Vt  =  -^(x0,t0),  Vxt  =  faft&O’to)’  XL  =  xRi  and  XR 

to  the  right  endpoint  of  the  next  cell  to  the  right,  where  Xq  refers  to  the  next  cell  to 
the  right,  whose  left  endpoint  is  x0;  tf  is  unchanged.  Go  to  step  23. 

17.  Like  step  16,  replacing  right, ,iCi,left  with  left, xl, Xq, x r, right,  respectively. 

18.  [tracking  reaches  right  endpoint  at  end  of  time  step]  Set  =  tf ,  x0  =  x(t0)  =  £'#, 

Ov 

v0  =  v(x0,  t0),  vt  =  —  (^o5^o )?  and  tf  to  the  end  of  the  next  time  step. 

\J  b 

Ov  O^v 

19.  If  u0  >  0,  or  if  u0  =  0  and  vt  >  0,  set  vx  =  —(x£,t0),  vxt  =  g^(xo^o),  xL  =  xR, 

and  xr  to  the  right  endpoint  of  the  next  cell  to  the  right.  If  w0  =  0  and  vt  <  0,  set 
Ov  O^v 

vx  =  —  (xq  ,  t0)  and  vxt  =  (xU  ,tT);  xl  and  xr  are  unchanged.  In  either  case,  go 

Ox  OxOt 

to  step  23. 

20.  Like  step  18,  replacing  right, xr  with  left,.x'L. 

21.  Like  step  19,  replacing  >,>,Xq  ,xl,xr, right, <,  Xq  with  <,<,Xq  ,xr,xl, left,>,o:o~. 

22.  [track  through  the  first  of  two  reversals,  staying  in  the  current  cell]  Set  to  —  t, 
xo  -  a; (t0),  U)  -  v(x0,to),  vx  -  —(x0,t0),  and  vt  -  ^(x0,t0);  vxt,  tf,  xL,  and  xR 


Ox 

are  unchanged.  Go  to  step  23. 

23.  [one  local  tracking  step  completed]  Go  to  step  1. 


4  RESULTS  AND  DISCUSSION 
4.1  Example  1 

We  first  consider  the  one-dimensional  problem  presented  in  detail  by  Lu  (1994).  The  sim¬ 
ulated  region  is  shown  in  Fig.  4  and  consists  of  an  aquifer  of  thickness  10  m  that  extends 
from  a  fully  penetrating  ditch  at  x  —  0  to  oo.  Aquifer  hydraulic  conductivity  is  .0002  m/s, 
specific  storage  is  2.0E-8  m-1,  and  porosity  is  0.5.  Initially  all  heads  in  the  aquifer  and  ditch 
are  equal  to  10  m.  At  time  t  —  0,  the  head  in  the  ditch  is  instantaneously  lowered  to  0  and 
maintained  at  that  level  for  all  subsequent  times.  Velocities  at  any  point  in  space  and  time 
can  be  determined  from  the  analytical  solution  given  by  Carslaw  &  Jaeger  (1959). 
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Figure  4:  Aquifer,  10  m  thick,  simulated  in  example  1  (Lu  1994). 

Table  1:  Particle  travel  times  (days)  for  example  problem  1  with  number  of  time  steps  NT 
=  15,  injection  time  T  (minutes)  =  1  or  1000. 


Method 

Travel  time 
T=1  T=1000 

Euler 

14.35 

21.15 

present 

13.28 

20.78 

Lu 

13.28 

20.77 

Pollock 

9.28 

17.43 

We  used  the  same  discretizations  as  Lu  (1994):  node  spacings  of  0.5  m  and  15  time 
intervals  (0.00035,  0.001,  0.01,  0.05,  0.2,  0.7,  1.2,  2.0,  3.0,  5.0,  9.0,  13.0,  17.0,  21.0,  30.0 
days).  Particles  were  tracked  to  the  ditch  from  a  distance  of  5  m  away  for  two  initial  times 
of  release  (1  and  1000  min). 

Travel  times  between  the  injection  point  and  the  ditch  are  presented  in  Table  1  (NT  =  15) 
for  Euler  integration  with  1000  steps  (Lu  1994),  the  method  proposed  here,  the  method  of 
Lu  (1994),  and  the  semianalytical  method  of  Pollock  (1989)  which  assumed  that  within  each 
time  interval  velocity  was  constant  and  equal  to  that  at  the  end  of  the  time  interval.  Values 
calculated  for  the  latter  two  methods  match  those  presented  by  Lu  (1994)  and  thus  imply 
that  our  computer  program  is  correctly  calculating  travel  times.  The  difference  in  results 
between  our  method  and  Lu’s  is  trivial  for  this  example  and  indicates  that  our  method  has 
no  advantage.  Indeed,  the  increased  computational  effort  required  by  our  method  makes  it 
less  attractive  than  Lu’s  for  this  particular  problem. 

In  order  to  test  sensitivity  of  results  to  time  discretization,  the  problem  was  rerun  with 
100  variably  spaced  time  steps  (NT  =  100).  Again  we  saw  virtually  no  difference  between  our 
method  and  Lu’s;  however,  results  of  all  methods  were  closer  to  the  Euler  solution  with  the 
finer  time  discretization.  In  particular,  the  results  of  Pollock  (1989)  appeared  more  attractive, 
especially  in  light  of  the  lesser  computational  requirements  it  has.  A  similar  sensitivity  was 
conducted  on  the  spatial  discretization;  it  showed  very  little  change  in  results  for  refined 
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grids. 

4.2  Example  2 

For  the  second  test  problem,  we  modify  the  first  example  by  extending  to  two  dimensions 
and  by  allowing  heterogeneity  in  the  hydraulic  conductivity  field.  The  simulated  domain  now 
extends  5  m  in  the  y-direction  with  spacing  equal  to  that  in  the  ^-direction  (0.5  m).  The 
ditch  extends  across  the  entire  length  of  the  grid  in  the  y-direction  and  initial  and  boundary 
conditions  remain  unchanged  so  that  flow  is  still  predominantly  in  the  ^-direction. 

The  natural  log  of  hydraulic  conductivity  was  randomly  generated  by  the  method  of 
turning  bands  (Mantoglou  &  Wilson  1982).  Spatial  correlation  was  assumed  over  a  distance 
of  3  m  and  was  described  by  a  spherical  variogram.  The  mean  and  variance  of  the  generated 
log  values  were  -3.912  and  5,  respectively.  Specific  storage  was  again  set  at  2.0E-8  m-1.  The 
flow  equation  was  solved  numerically  with  the  grid  extended  in  the  ^-direction  to  a  distance 
sufficient  to  insure  no  boundary  effects.  The  boundaries  at  y  =  0  and  5  m  were  impermeable. 

At  an  initial  time  of  200  sec,  100  particles  were  introduced  at  x  =  5  m  (10  particles  equally 
spaced  in  each  of  the  10  cells  in  that  column  of  the  grid).  Particles  were  tracked  through  the 
complex,  2-dimensional  flow  field  until  they  reached  the  ditch.  For  each  particle  the  starting 
and  ending  location  and  the  total  travel  time  were  recorded.  Results  of  different  methods 
were  compared  on  the  basis  of  travel  time  and  ending  location.  For  each  realization  of  the 
log  hydraulic  conductivity  field,  three  simulations  were  made.  These  simulations  differed 
only  in  the  time  discretization  (number  of  time  steps  NT  =  200,  100,  or  20).  Results  from 
many  simulations  were  studied  and  found  to  be  qualitatively  very  similar.  Therefore  we 
present  results  from  only  two.  To  analyze  results  we  consider  the  fine  time  discretization 
(NT  =  200)  using  the  present  method  to  be  the  standard  for  comparison.  This  decision 
is  somewhat  justified  by  the  fact  that  there  was  virtually  no  difference  in  results  among 
the  three  methods  for  this  discretization,  and  the  present  method  accounts  for  the  highest- 
order  velocity  variations.  For  each  particle  we  compared  the  total  travel  time,  t?,  and  final  y 
coordinate,  yf,  to  the  values  obtained  from  the  standard  solution  (ts.  ys,  respectively).  The 
mean  square  error  (MSE)  of  travel  time  was  calculated  from  normalized  values  of  (ts—t^)/ts. 
The  average  time  deviation  (DEV)  and  y-location  deviation  were  calculated  as  average  values 
of  \ts  —  tJ | / ts  and  | ys  —  y*\/ys.  Values  for  these  terms  are  contained  in  Table  2. 

The  significance  of  the  enhanced  accuracy  of  the  present  method  can  be  measured  in 
two  ways:  by  the  differences  between  its  results  and  those  the  other  methods,  and,  when 
NT  is  decreased,  by  the  amounts  by  which  the  errors  of  the  other  methods  exceed  those 
of  the  present  method.  For  the  fine  time  discretization  (NT  =  200),  there  is  very  little 
difference  between  results  of  our  method  and  that  of  Lu  (1994),  while  the  Pollock  (1989) 
method  showed  somewhat  larger  differences.  As  NT  was  decreased,  the  differences  among  all 
methods  increased,  and  this  trend  toward  differences  of  the  order  of  1%  or  more  was  apparent 
in  all  of  the  realizations  examined.  When  NT  changes,  so  does  the  computed  velocity  field; 
this  causes  the  largest  share  of  the  errors  in  the  bottom  rows  of  the  table.  However,  the 
present  method  appears  to  reduce  the  NT-sensitivity  of  the  method  of  Lu  (1994)  by  the 
order  of  10%  for  this  problem. 
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Table  2:  Results  of  example  2  for  number  of  time  steps,  NT,  equal  to  200,  100,  20. 


Realization  1 

Realization  2 

Method 

Average 
travel  time 

Time 

MSE 

Time 

DEV 

y 

DEV 

Average 
travel  time 

Time 

MSE 

Time 

DEV 

y 

DEV 

NT  =  200 
present 

8.98E+5 

_ 

_ 

_ 

1.64E+5 

_ 

_ 

_ 

Lu 

8.98E+5 

0.000 

0.028 

0.001 

1.64E+5 

0.000 

0.002 

0.001 

Pollock 

8.88E+5 

0.002 

0.036 

0.000 

1.60E+5 

0.001 

0.033 

0.001 

NT  =  100 
present 

8.93E+5 

0.007 

0.067 

0.001 

1.35E+5 

0.013 

0.090 

0.001 

Lu 

8.92E+5 

0.007 

0.071 

0.001 

1.35E+5 

0.014 

0.094 

0.002 

Pollock 

8.79E+5 

0.015 

0.120 

0.001 

1.28E+5 

0.027 

0.153 

0.001 

NT  =  20 
present 

9.24E+5 

0.064 

0.249 

0.001 

1.37E+5 

0.028 

0.165 

0.002 

Lu 

9.09E+5 

0.081 

0.271 

0.005 

1.35E+5 

0.032 

0.177 

0.005 

Pollock 

8.84E+5 

0.110 

0.325 

0.001 

1.26E+5 

0.056 

0.237 

0.001 

Table  3:  Results  of  example  3  for  number  of  time  steps,  NT,  equal  to  200,  100,  20. 


Realization  1 

Realization  2 

Method 

Average 
travel  time 

Time 

MSE 

Time 

DEV 

y 

DEV 

Average 
travel  time 

Time 

MSE 

Time 

DEV 

y 

DEV 

NT  =  200 
present 

2.87E+6 

_ 

_ 

_ 

1.20E+7 

_ 

_ 

_ 

Lu 

2.93E+6 

0.003 

0.019 

0.000 

1.20E+7 

0.000 

0.005 

0.001 

Pollock 

2.79E+6 

0.001 

0.034 

0.000 

1.18E+7 

0.009 

0.047 

0.001 

NT  =  100 
present 

2.60E+6 

0.006 

0.053 

0.000 

1.01E+7 

0.013 

0.080 

0.000 

Lu 

2.62E+6 

0.009 

0.061 

0.008 

1.00E+7 

0.014 

0.081 

0.002 

Pollock 

2.57E+6 

0.008 

0.069 

0.007 

9.93E+6 

0.033 

0.131 

0.000 

NT  =  20 
present 

2.87E+6 

0.040 

0.164 

0.007 

1.00E+7 

0.024 

0.134 

0.002 

Lu 

2.14E+6 

0.048 

0.177 

0.015 

9.83E+6 

0.032 

0.147 

0.005 

Pollock 

2.21E+6 

0.061 

0.201 

0.020 

1.02E+7 

0.068 

0.223 

0.013 

4.3  Example  3 

The  third  test  problem  is  identical  to  the  second,  except  that  ground- water  velocities  change 
direction  during  the  simulation.  Water  level  in  the  ditch  is  initially  lowered  to  0  m,  then 
raised  back  up  to  10  m,  and  then  lowered  to  0  m  again.  This  causes  the  prominent  direction 
of  flow  to  be  first  toward,  then  away,  and  finally  back  toward  the  ditch.  This  example  should 
produce  relatively  large  values  for  d2v/dxdt  and  should  therefore  be  indicative  of  the  largest 
differences  to  be  expected  between  our  method  and  that  of  Lu  (1994).  Results  of  these 
runs  are  shown  in  Table  3.  The  reduction  of  the  NT-sensitivity  is  somewhat  greater  than  in 
Example  2,  and  the  differences  in  results  in  general,  and  particularly  for  the  coarsest  time 
discretization  in  Realization  1,  are  notably  greater.  In  this  regime,  inclusion  of  the  bilinear 
term  can  lead  to  highly  significant  differences  in  particle  travel  times. 
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5  CONCLUSIONS 

The  new  particle-tracking  method  produces  analytical  trajectories  in  velocity  fields  that 
arise  from  a  lowest-order  Raviart-Thomas  spatial  discretization  on  a  cartesian  grid,  fully 
accounting  for  linear  temporal  variations,  including  space-time  bilinearity.  Inclusion  of  the 
bilinearity  is  unimportant  in  some  cases,  but  with  spatial  heterogeneity,  and  particularly 
with  time- varying  pumping  or  recharge,  it  can  make  large  differences  in  travel  times. 
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